Assessing the rainfall infiltration on FOS via a new NSRM for a case study at high rock slope stability

Assessing the stability characteristics of high rock slope under rainfall via theoretical research, numerical simulation, and field monitoring is of great implications for safety construction in open-pit mine engineering. Thus, based on the Hoke-Brown criterion, instantaneous internal friction angle and cohesion of high-slope rock mass under high stress conditions were deduced, and a nonlinear strength reduction method for high rock slope was established. The safety factors of the open-pit mine were calculated by COMSOL Multiphysics, which considering the high rock southwest slope and detected rainfall in Dagushan Open-pit Mine, China. The results showed that high rock slope stability could be more accurately analyzed by the proposed method due to its full consideration of slope stress state effect compared with the equivalent Mohr- Coulomb method. When the slope is low, the difference between the calculation results of the equivalent Mohr- Coulomb method and the proposed method is small, but with the increase of the slope height, the difference between the two calculation results gradually increases. When the transient saturated is formed in the slope surface layer and gradually increases, the reduction rate of the factor of safety (FOS) gradually increases. When the total rainfall is the same, the effect of short-term heavy rainfall on slope stability is less than that of long-term ordinary rainfall. The results obtained form this work provided important insights into the stability of high rock slope and rainfall infiltration in open-pit mine, and the safety factor is crucial for guiding the mining process design.

In recent decades, with the rapid development of computer technology, numerical analysis method has been widely used in slope engineering [13][14][15][16] . On this basis, the strength reduction method has gradually become the main method to study the stability of rock slope. Compared with other methods, the rock stress-strain relationship could be considered by this method and that obtain the critical sliding surface and the minimum safety factor without assuming the shape and position of the sliding surface [17][18][19][20] . In the meantime, the progressive failure process of slope could be simulated by the strength reduction method. The generalized Hoek-Brown has been proved to be effective in the empirical calculation of shear strength of rock mass 21,22 . Therefore, it is necessary to combine the strength reduction method with Hoek-Brown criterion [23][24][25] . Many researchers have worked on the study and implementation of strength reduction method (SRM) with the nonlinear generalized Hoek-Brown (GHB) criterion, referred to as the nonlinear strength reduction method (NSRM) [26][27][28][29] , as shown in Table 1. Hammah et al. 30 proposed a global GHB shear strength reduction method by lowering the GHB shear envelopes. Benz et al. 31 established a mathematical expression for the relationship between the GHB criterion and the MC criterion, and used the MC criterion reduction factors to determine the strength reduction based on GHB criterion. Wu et al. 32 concluded that for intact rock, only uniaxial compressive strength and rock hardness should be used as strength indicators for GHB strength reduction. Fu et al. 33 introduced the Mohr-Coulomb instantaneous friction angle as a variable, used the instantaneous Moore-Coulomb friction angle and cohesive strength to describe the shear strength of rock mass, and proposed a nonlinear shear strength reduction technique for slope stability calculations. Sun et al. 34 calculated the factor of safety (FOS) by combining slope geometric parameters, rock mass property parameters and GHB parameters, and proposed a simplified method to realize nonlinear strength reduction based on GHB criterion.
In summary, the nonlinear strength reduction methods have been proposed by previous researchers from different perspectives, but a strong focus on the slope stress state is missing. In practical engineering, especially in the application of high rock slope, the stress variation range is large, and the strength parameters of rock mass need to be reasonably determined for strength reduction. For high rock slope, the high stress state of rock mass is the most critical factor for strength reduction. Therefore, a nonlinear strength reduction method considering high stress state should be proposed for the FOS calculation. These knowledge gaps are seen as the main objectives in this work of stability analysis on the high rock slope. Aiming at the high stress of high rock slope, the instantaneous internal friction angle and instantaneous cohesion of rock mass under different stress states are deduced, and the new nonlinear strength reduction method for high rock slope is established according to the relationship between normal stress and shear stress of rock mass under the Hoke-Brown criterion. Taking the southwest slope of Dagushan Iron Mine as the research background, the safety factors of high rock slope under different rainfall conditions are calculated by the proposed nonlinear strength reduction method on COMSOL Multiphysics. The influence of rainfall infiltration on the stability of high rock slope in open-pit mine is analyzed.

A new nonlinear strength reduction method
Hoek-Brown yield criterion. Hoek-Brown (HB) strength criterion is an empirical formula of rock failure based on the theoretical research results of rock mechanics by evert Hoek and E.T. Brown. They have analyzed a large number of engineering field data. The theory scientifically sums up a large number of in-situ test results, solves many practical engineering problems, and is widely used in open-pit mining and slope engineering.
Hoek et al. 14 introduced the geological strength index GSI and disturbance coefficient D of the rock mass into the theory. The modified generalized HB expression is as follows: where m i is Hoek Brown constant of complete rock; GSI is geological strength index of rock mass, which is selected according to rock mass structure and surface characteristics of structural plane; D is disturbance coefficient, which is 1.0 for conventional blasting, 0.7 for mechanical excavation and 0.8-0.9 for controlled blasting.
Equivalent Mohr-Coulomb strength parameters. Because most of the geotechnical engineering software is still based on Mohor-Coulomb (MC) failure criterion, few has been studied on non-linear shear strength reduction techniques [35][36][37][38][39][40] . It is necessary to determine the cohesion and internal friction angle of rock mass by HB criterion. HB criterion is non-linear, while MC criterion is linear, so cohesion and internal friction angle cannot be calculated directly [41][42][43][44][45] . Therefore, the strength curve of HB criterion is fitted as equivalent MC straight line to obtain cohesion and internal friction angle of rock mass, as shown in Fig. 1.
In the fitting process, the upper and lower difference areas of MC straight line are balanced to obtain the equivalent MC parameters.
According to the actual situation of slope engineering, Hoke et al. made specific provisions for the parameters in Formula (5)- (7): where σ cm is the compressive strength of rock mass, σ 3max is the maximum confining pressure, H is the height of slope and γ is the weight of rock mass. As shown in Fig. 2, in the coordinate space with normal stress σ as the horizontal axis and shear stress τ as the vertical axis, the equivalent MC strength envelope curve is a straight line and HB strength envelope is a kinked line. The instantaneous cohesion c i and internal friction angle φ i can be obtained by making a tangent line on the HB strength envelope curve 33,35 . It can be seen from the Fig. 2 that with the increase of stress, the cohesion c i increases gradually, then it will be equal to the equivalent cohesion at a certain stress, while the internal friction angle φ i decreases gradually, then it will also be equal to the internal friction.
In slope engineering, the stress is relatively small in the shallow part of the slope as shown in point 1 in Fig. 2, the cohesion and internal friction angle can be C 1 and φ 1 by tangent line, C 1 is smaller than C, and φ 1 is larger than φ. However, in the deep part of the slope, the stress is relatively large as shown in point 2 in Fig. 2, the cohesion and internal friction angle can be obtained by tangent line, which are C 2 and φ 2 respectively, C 2 is larger than C, and φ 2 is less than φ.
Therefore, in view of the high stress characteristics of high rock slope, the establishment of the relationship between cohesion, internal friction angle and minimum principal stress is the key of slope strength reduction calculation.
Nonlinear strength reduction method. The equivalent MC strength parameters are fitted in the range of σ t < σ 3 < σ 3max . When the minimum principal stress is small, the calculated strength parameters are relatively accurate, but the larger the minimum principal stress is, the lower the calculation accuracy is. Because of the high and steep rock slope in open pit, it is necessary to determine the strength parameters of rock mass in case of large stress.
Under the condition of complex stress, the stress state of yield point in rock mass can be represented by limit Mohr circle on the coordinate space with normal stress σ as horizontal axis and shear stress τ as vertical axis. Whether the strength envelope is a straight or kinked line, only one point of the limit Mohr circle is tangent to the envelope curve, as shown in Fig. 3.
The cohesion and internal friction angle corresponding to any point on the HB strength envelope curve are expressed by c i and φ i respectively. The normal stress σ and shear stress τ on this point can be expressed by MC theory as follows:   For HB failure criterion, because the failure envelope curve is a kinked line, φ i changes with the change of minimum principal stress state, that is, Taking the derivative of Formula (12) and Formula (13) with respect to σ 3 yields: On the τ-σ coordinate space: Substituting Formula (14) and Formula (15) into Formula (16) and simplifying yields: The equation of MC failure envelope curve on τ-σ coordinate space can be expressed as follows: Substituting Formula (12), (13) and (17) into Formula (18) yields: It can be seen from Formula (17) and (19) that the values of c i and φ i are related to the minimum principal stress in the limit state, the HB characteristic parameters of rock mass and the uniaxial compressive strength of rock. The HB characteristic parameters and uniaxial compressive strength of HB rock mass are known. The stress state of the slope rock mass is calculated by the finite element method, and the minimum principal stress σ e3 is obtained. As shown in Fig. 4, c i and φ i of each point in the rock mass can be obtained by substituting σ 3 = σ e3 .

Implementation of nonlinear strength reduction method based on HB criterion. Using COM-
SOL multiphysics software, this paper embeds the calculation method of instantaneous equivalent MC strength parameters into the strength reduction calculation program, so as to realize the application of HB criterion finite element strength reduction method in slope stability analysis. The calculation flow is shown in Fig. 5. The calculation process consists of three parts: establishing models, parameters calculation and strength reduction. In the calculations, the instantaneous internal friction angle and instantaneous cohesion are derived according to Hoek-Brown criterion, while the elastic-plastic finite element calculation of slope and minimum principal stress are calculated according to Mohr-Coulomb criterion. The proposed method includes the following steps: Step 1 The solution area is determined according to the geometric parameters of slope, and the material parameters of rock mass are inputted. Mesh generation is carried out according to certain laws in the solution area, and the triangular element with three nodes is selected. The loading and boundary conditions are set, and the initial slope stability can be calculated.
Step 2 The initial value of reduction coefficient K is set to 1, the stress state of slope rock mass is calculated by FEM. The minimum principal stresses for each element are obtained as the initial conditions for the elastoplastic calculations. Substituting the known Hoek-Brown parameters and the minimum principal stress into  (17) and (19), the instantaneous cohesion ci and internal friction angle φi of each element of the slope model are calculated.
Step 3 According to the reduction coefficient K, ci and φi are reduced. Judge whether the slope model converges. If the model converges, end the calculation and output the reduction coefficient K. Otherwise, according to the reducted parameters, the Mohr-Coulomb criterion is used for the elastic-plastic finite element calculation of slope to obtain the minimum principal stress of each element, and then instantaneous strength

Case study and verification
In order to verify the effectiveness of the new nonlinear strength reduction method proposed in the previous section, a simple slope model is established using COMSOL multiphysics software, with a slope height of 200 m and a slope gradient of 1:1, as shown in Fig. 6. The model boundary conditions are set as followed: the left and right boundaries of the slope are set as roller support constraints, and the bottom boundary is set as fixed constraints. Case A adopts the new nonlinear strength reduction method proposed in the previous section, and case B adopts the equivalent MC strength parameter reduction. The strength parameters of slope rock mass are based on the data in Table 2, and the strength parameters of equivalent MC criterion are calculated according to Formula (5) and Formula (6). The stress state in the slope is calculated using COMSOL multiphysics software, and the calculation result of the minimum principal stress σ 3 is shown in Fig. 7.
The calculated minimum principal stress σ 3 and HB strength parameters are substituted into Formula (17) and Formula (19) to obtain the non-linear MC strength parameters c i and φ i in the slope. The calculation results are shown in Fig. 8.   Fig. 8, it can be seen that the cohesion of case A decreases gradually from the inside out, and the cohesion at the upper part is less than the equivalent cohesion, while the change of internal friction angle is opposite to the cohesion, they are non-linear.
For the case A, the strength reduction calculation is carried out by using the nonlinear MC strength parameters c i and φ i based on HB criterion. The safety factor of the slope is 1.24. For the case B, the equivalent MC criterion strength parameters are directly used for strength reduction calculation. The safety factor of the slope is 1.30. The plastic strain distribution and total displacement diagram of the two cases in the limit state are shown in Fig. 9.
It can be clearly observed from Fig. 9 that the potential sliding surface of case A is closer to the slope surface than case B, and the landslide body of case A is smaller than case B. And the factor of safety (FOS) calculated in case A is less than that calculated in case B and the stability of slope in case B is overestimated. Obviously, the simulation results of case A are more closer to the actual situation. Therefore, the proposed nonlinear strength reduction method is more reasonable in numerical calculation.
The equivalent MC strength parameters are estimated by fitting, and their accuracy depends on the range of the minimum principal stress. However, the nonlinear strength parameters derived in the previous section are related to the minimum principal stress, the characteristic parameters of HB rock mass and the uniaxial   In order to further explore the rationality and practicability of the proposed nonlinear strength reduction method for the stability calculation of high rock slope, the authors established rock slope with different height in the COMSOL multiphysics software, and analyzed the difference of the factor of safety (FOS) calculated using the two methods with the change of slope height, as shown in Fig. 10. Figure 10 shows the relationship between safety factor and slope height. It can be seen from the figure that as the height increases, the safety factor of the whole slope decreases gradually. In the meantime, when the slope is low, the difference between the two safety factors is small, but with the increase of the slope height, the difference between the calculation results of the equivalent Mohr-Coulomb method and the nonlinear strength reduction method gradually increases. The characteristics of high rock slope with a wide range of stress variation is not fully considered in the equivalent Mohr-Coulomb method and the slope stability calculation results may be overestimated. Therefore, the nonlinear strength reduction method based on HB criterion could be effectively applied to the stability analysis of high rock slope.

Application on southwest slope of Dagushan Open-pit Iron Mine
The effect mechanism of rainfall infiltration. Rainfall infiltration is the process of rainwater infiltration through the surface into rock mass, which is essentially the movement process of rainwater in the porous medium constantly displacing air. The saturated zone is below the water table, and the shallow surface layer will evolve into a transient saturated zone under the influence of rainfall. The unsaturated zone is between transient saturated zone and water table, as shown in Fig. 11. Overall, the effect mechanism of rainfall infiltration on the stability of rock slope includes two aspects. On the one hand, rainfall can cause the groundwater level to rise, leading to a transient saturated zone in the surface layer of slope and raising the pore-water pressure in the corresponding area. On the other hand, rock mass can be affected by the softening effect of rainwater. Rainfall

The establishment of geometric model. Based on the field investigation of Dagushan Open-pit Iron
Mine and the collection of relevant geological data and monitoring data, the southwest slope area is selected as the research object, which is a typical high rock slope. The main rock layers of the slope in this area are migmatite, porphyrite, granite and magnetite quartzite. The distribution of the rock layers in the section is shown in Fig. 12. The stratum in this area is composed of many kinds of rock strata, and the joint fissures are relatively developed. Through communication with the site staff, there has been a small-scale slide on the slope in the area, which causes the deformation of the upper part of the slope and damages the corresponding facilities and buildings. According to the microseismic monitoring data, a large number of microseismic events occurred in this area due to rainfall. In this work, COMSOL multiphysics multi physical field simulation software is used for numerical calculation, and Richards equation model is used for seepage field calculation. The selected calculation profile of the southwest slope of Dagushan open-pit iron mine is shown in Fig. 12. The whole calculation grid area is 411,500 square meters, the unit type is triangle unit, the number of units is 3245, and the number of grid vertices is 1719.

Initial conditions and parameters. Permeability tensor and hydraulic characteristic parameters of frac-
tured rock mass. Through the field sampling and investigation of the joints and fissures in the slope rock mass of Dagushan open-pit mine, and then statistical analysis, the fracture geometric parameters of the southwest slope rock mass are obtained as shown in Table 3.   Table 2, the equivalent permeability tensor of the corresponding rock mass is calculated by the fracture sampling measurement method, and the calculation results of the equivalent permeability tensor of each rock layer are shown in Table 4.
The water retention curves are fitted by the Van Genuchten model, which can be expressed as: where is the saturation; θ w is the moisture rate of rock mass; θ s is the saturated moisture rate of rock mass; θ r is the residual moisture rate of rock mass; ψ is the matrix suction; α , n and m are characteristic parameters, and m = 1 − 1 n. The fitted hydraulic characteristics of rock slope are shown in Table 5.
Mechanical parameters of rock mass. Through field investigation and data collection, the development of joint fissures, dominant joint occurrence and structural plane distribution characteristics of slope rock mass can be obtained. According to the rock mechanics indexes obtained from various rock mechanics tests conducted before in Dagushan, the calculation parameters of slope stability are obtained based on HB strength theory. The adopted values of relevant calculation parameters are shown in Table 6.
Rainfall schemes. The change of rainfall conditions will affect the distribution of seepage field in the rock slope, and the seepage field will affect the stability of the slope. In this work, we will study the change of slope seepage field and then calculate the change of slope stability under different rainfall conditions. According to the relevant meteorological data, the daily maximum precipitation of Dagushan open-pit mining area is 236.8 mm.
On this basis, to analyze the influence of rainfall intensity, rainfall duration and total rainfall on slope stability, we designed four different schemes as shown in Table 7.   Table 6. Parameters for stability calculation.

Uniaxial compression stress (MPa) Elastic modulus (GPa) Poisson's ratio
Hoek www.nature.com/scientificreports/ Initial conditions. The boundary of the left and right sides of the slope numerical model is set with a given water head boundary, the water head is equal to the height of the groundwater level at the boundary. In the meantime, the bottom boundary is set as flux boundary due to hydraulic conduction. The boundary conditions of the slope surface are set according to the rainfall intensity: when the rainfall intensity is less than the saturated permeability coefficient of the slope rock, the surface boundary is set as a given flow boundary, and the flow size is rainfall intensity; if not, the boundary is set as the given water head boundary, and the surface runoff formed at this time is relatively small, so the water head height is the elevation of this place, that is, the surface of the slope is the zero-pressure surface. The initial conditions are obtained from the hydrogeological conditions of the region according to the back analysis of the stable seepage field, and the pore water pressure distribution in line with the actual situation is obtained by continuously adjusting the flow flux at the bottom boundary. The calculated initial pore water pressure distribution is shown in Fig. 13, in which the blue solid line is the groundwater level line, the pore water pressure in the area below the groundwater level line is positive, and the pore water pressure in the area below the groundwater level line is negative, and the negative value represents the matrix suction in the area.
Seepage simulation analysis of high rock slope under rainfall infiltration. The four rainfall schemes above are numerically simulated by COMSOL Multiphysics. Within the solution area, 1-1 characteristic profile, 2-2 characteristic profile and 3-3 characteristic profile are selected at the top, middle and bottom of high rock slope for analysis, and the locations of characteristic profiles are shown in Fig. 14. And the calculation results are shown in Fig. 15, 16, 17, 18, 19, 20, 21, 22. The simulation results of the four rainfall schemes show that: (1) Under continuous rainfall, unsaturated zone on the surface layer of slope gradually evolves into a transient saturated zone and a positive pore-water pressure Total rainfall (mm) 120 240 240 240 Figure 13. Initial pore water pressure and groundwater level. www.nature.com/scientificreports/ appears. When the rainfall intensity is constant, the longer the rainfall time is, the larger the transient saturation zone and the hydraulic conductivity depth is. (2) During rainfall, the pore-water pressure of the surface layer of slope gradually increases, and the wetting front evolves to the deeper part. After the rainfall stops, the pore-water pressure gradually decreases in the area above a certain depth; in the area below this depth, the pore-water pressure gradually increases, and the wetting front continues to evolve to a deeper part. (3) When the total rainfall is constant, the greater the rainfall intensity, the greater the growth rate of pore-water pressure in the surface layer of slope, and the greater the maximum pore-water pressure. The longer the rainfall duration, the greater the hydraulic conduction depth. After the rainfall stops, the increase of hydraulic conduction depth is basically the same. (4) During the calculation time, the rainfall infiltration has little effect on the groundwater level because the hydraulic conductivity has not reached the water table surface of groundwater.

Stability analysis of high rock slope under rainfall infiltration
The proposed nonlinear strength reduction method based on Hoek Brown criterion is used to calculate the safety factor of slope. The minimum principal stress in the slope is calculated so that the cohesion and internal friction angle of the slope rock mass can be obtained. The pore water pressure of slope at the initial time of rainfall (0 h) is substituted into the calculation model, and the strength reduction calculation is carried out for the cohesion and internal friction angle calculated above. Based on the non convergence of calculation, the factor of safety (FOS) at the initial time of rainfall is calculated to be 1.544. The distribution of plastic strain and total displacement at the initial time are shown in Fig. 23.  www.nature.com/scientificreports/ Influence of rainfall duration on slope stability. The rainfall intensity of Scheme I and Scheme II is 10 mm/h. The rainfall duration of Scheme I is 12 h, the rainfall duration of Scheme II is 24 h, and the calculated duration after the rainfall stops is 24 h. According to the calculation results of the slope safety factor of Scheme I and Scheme II, the change diagram of the slope safety factor with time can be obtained, as shown in Fig. 24. It can be seen from Fig. 24 that the safety factor of the slope is 1.544 at the beginning of the rainfall. When the rainfall intensity remains unchanged, the safety factor decreases gradually with the increase of rainfall duration. The safety factor is 1.443 and decreases by 0.101 in 12 h of rainfall. The rainfall of Scheme I stops at this time. Then the safety factor starts to rise gradually, and it rises to 1.504 at 24 h after stopping the rainfall. This is because the increment of water load is large, the pore water pressure is increasing and the matrix suction is decreasing when it rains. It leads the safety factor of slope to decrease. When the rainfall stops, there is no rainwater supply, and  www.nature.com/scientificreports/ the previously infiltrated rainwater continues to infiltrate under the effect of gravity. Then the pore water pressure on the upper part of the slope surface gradually decreases. It leads the safety factor to increase. But because the dissipation speed of pore water pressure is relatively slow, the safety factor gradually tends to be stable but less than the safety factor at the initial time.
In Scheme II, the rainfall continues for 24 h and the safety factor decreases to 1.412 at 24 h after raining. After 16 h, the decrease rate of safety factor increases, which is due to the formation of transient saturated zone (TSZ) on the surface of the slope. TSZ has a great impact on the stability of the slope. After the rainfall stops, the safety factor of Scheme II rises gradually to 1.483 at 24 h which is smaller than the safety factor of Scheme I at 24 h after rain.  (1) The same rainfall duration The rainfall intensity of Scheme I is 10 mm/h, and that of Scheme II is 20 mm/h. The rainfall duration of the two schemes is the same, both of them are 12 h, and the calculated time length after the rainfall stops is 24 h. According to the calculation results of the slope safety factor of Scheme I and Scheme III, the change diagram of the slope safety factor with time can be obtained, as shown in Fig. 25.
It can be seen from Fig. 25 that the changes in safety factor of Scheme I and Scheme III are similar, both of them reach the lowest level in 12 h of continuous rainfall, and the safety factor of Scheme III is lower than that Scheme I. In the meantime, in the duration of rainfall and the time after the rainfall stopped, the safety factor of Scheme III is always lower than that of Scheme I. For Scheme III, the decrease rate of safety factor increases   www.nature.com/scientificreports/ after 8 h, mainly because the formation of TSZ has a greater impact on slope stability. Overall, the comparison of Scheme I and Scheme III shows that when the rainfall duration is the same, the greater the rainfall intensity is, the lower the slope safety factor is (2) The same total rainfall The total rainfall of Scheme II, Scheme III and Scheme IV is 240 mm, the rainfall intensity is respectively 10 mm/h, 20 mm/h and 40 mm/h, the duration is respectively 24 h, 12 h and 6 h, and the calculated duration after the rainfall stops is 24 h. According to the calculation results of the slope safety factor of Scheme II, Scheme III and Scheme IV, the change diagram of the slope safety factor with time can be obtained, as shown in Fig. 26.
It can be seen from Fig. 26 that the safety factors of the three schemes gradually decrease with the continuous rainfall. All of the safety factors become the smallest at the end of the rainfall, and then they gradually recover and increase in a period of time after the rainfall stops. In the process of continuous rainfall, the safety factor of Scheme IV decreases the fastest, followed by Scheme III. This is due to the large rainfall in a short time, the large increment of water load in the slope and the large change of safety factor. Although the total rainfall of the three schemes is equal, the minimum safety factors of the three schemes are different. The minimum safety factor of Scheme IV is the largest, Scheme III is the second, and Scheme II is the smallest. When the rainfall stops (24 h), the safety factor of Scheme 4 is also the largest among the three schemes. The main reason for this phenomenon is that for heavy rainfall in a short time, although the water load increment is large, due to the short infiltration time of rainwater into the slope, the shallow hydraulic conduction depth, and some rainwater will form surface runoff. Although the rainfall intensity for a long time is smaller, the rainwater has sufficient time to continuously infiltrate, the infiltration depth into the slope is greater, the area of rainfall infiltration is larger, and the impact on slope stability is greater. The comprehensive analysis of the three schemes shows that when the total rainfall is the same, the impact of short-term heavy rainfall on slope stability is less than long-term ordinary rainfall, that is, the longer the rainfall duration, the greater the impact on slope stability.
At the end of rainfall, the plastic strain distribution of slope are shown in Fig. 27. In the figure, (a), (b), (c) and (d) represent the plastic strain of Scheme I at 12 h, Scheme II at 24 h, Scheme III at 12 h and Scheme IV at 6 h respectively. Comparing Fig. 23 with Fig. 27, it can be seen that the overall plastic zone does not change much, but there is a through plastic zone on some single steps of the slope, especially in Scheme II, which shows that the impact of rainfall infiltration on the stability of some single steps is greater than that of the overall slope. In the process of rainfall, the pore water pressure on the slope surface increases continuously, forming a transient saturation zone, and the rock mass of some steps reaches a transient saturation state. At this time, the water content of rock mass increases, the shear strength decreases, and the sliding force increases, resulting in a significant reduction in the stability of some single steps. This situation often occurs in practical engineering, and long-term rainfall could cause the deformation and failure of a single step or several adjacent steps of the slope.
As mentioned above, a comparative analysis of the four schemes reveals that the FOS will decrease with the increase of rainfall duration. Under the same rainfall duration, the higher the rainfall intensity, the lower the FOS of slope. And under the same total rainfall, the effect of short-term rainstorm on slope stability is smaller than that of long-term ordinary rainfall, and the longer the duration of rainfall, the greater the impact on slope stability.

Discussion
The high rock slope stability analysis of Dagushan Mine under rainfall infiltration is studied and simulated by the proposed nonlinear strength reduction method. By comparing with the equivalent Mohr-Coulomb method, it is found that high rock slope stability could be more accurately calculated by the proposed method due to its full consideration of slope stress state effect. The seepage characteristics and stability of high rock slope under rainfall infiltration are simulated, and it is found that the effect of short-term heavy rainfall on slope stability is less than that of long-term ordinary rainfall under the same rainfall amount. Model applicability to stability studiy of high rock slope in open-pit mine is confirmed, and its simulation in high rock slope should be applied www.nature.com/scientificreports/ more widely. In this paper, the equivalent continuous medium model is used as the seepage model. However, the actual fracture yield, spacing and size are randomly distributed. It is necessary to determine the permeability tensor according to the fracture distribution characteristics to make the calculation more realistic. The numerical model is selected as a two-dimensional characteristic profile, while the geometry and rock structures of slope in engineering are anisotropic. The main limitation of the model is that it is only applicable to two-dimensional conditions. It is necessary to optimize the proposed model so that it can be applied to the analysis under threedimensional conditions in further research.

Conclusion
(1) According to the characteristics of high rock slope with large stress variation range, the instantaneous cohesion and internal friction angle of each slope element under different minimum principal stresses are derived based on Hoek-Brown strength criterion by considering the influence of slope element stress state. On this basis, the nonlinear strength reduction method for high rock slope is proposed. (2) By comparing with the equivalent Mohr-Coulomb method, it is found that the proposed method calculates a smaller landslide body and the potential sliding surface of slope is closer to slope surface. When the slope is low, the difference between the calculation results of the equivalent Mohr-Coulomb method and the proposed method is small, but with the increase of the slope height, the difference between the two calculation results gradually increases. (3) Under continuous rainfall, the pore-water pressure gradually rises and the matrix suction gradually decreases, causing the FOS to gradually decrease. When the rainfall stops, the pore-water pressure in the upper of slope surface layer gradually decreases, and FOS gradually increases but is smaller than the initial FOS. The transient saturated zone has a greater influence on the slope stability. When the transient saturated is formed in the slope surface layer and gradually increases, the reduction rate of FOS gradually increases. (4) Keeping the rainfall intensity, the greater the rainfall duration, the smaller the FOS; keeping the rainfall duration, the greater the rainfall intensity, the smaller the FOS. When the total rainfall is the same, the effect of short-term heavy rainfall on slope stability is less than that of long-term ordinary rainfall. (5) The proposed model could be widely used in the stability analysis of high rock slope, transient saturated zone and long-term ordinary rainfall should be treated as important concerns in policy implications of slope stability analysis. However, the main limitation is that it is not applicable to three-dimensional condition and simplifies the calculation of infiltration tensor, more simulation studies will be presented in our further work.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.